Phase-Field Formulation for Quantitative Modeling of Alloy Solidification 
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A phase- field formulation is introduced to simulate quantitatively microstructural pattern formation 
in alloys. The thin-interface limit of this formulation yields a much less stringent restriction on the 
choice of interface thickness than previous formulations and permits to eliminate non-equilibrium 
effects at the interface. Dendrite growth simulations with vanishing solid diffusivity show that both 
the interface evolution and the solute profile in the solid are well resolved. 
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The phase-field approach has emerged as a method of 
choice to simulate microstructural evolution during solid- 
ification [0-1 . This approach has the well-known advan- 
tage that it avoids to explicitly track a sharp boundary by 
smearing the interface region over some thickness ~ W . 
However, it has the disadvantage that it is hard to use 
quantitatively. This is because it is often computation- 
ally too stringent to choose W small enough to resolve 
the desired sharp- interface limit of the phase-field model, 
even on computers of today. This is especially true for 
small growth rates where the scale of the microstructure 
is typically several orders of magnitude larger than the 
microscopic capillary length do- 

Progress has recently been made to overcome this dif- 
ficulty by using a "thin-interface" analysis of the phase- 
field model ||-||] where W is assumed small compared to 
the scale of the pattern, but not smaller than do- For the 
standard symmetric model (with equal thermal conduc- 
tivities in the solid and liquid) , Karma and Rappel (KR) 
have shown that this thin-intcrfacc limit yields two es- 
sential benefits Q . Firstly, it maps onto the standard set 
of sharp-interface equations that one obtains in the clas- 
sical sharp- interface limit where W/do 0, but yields a 
much less stringent restriction on W/ do that renders the 
computations tractable. Secondly, it makes it possible to 
eliminate interface kinetic effects by a specific choice of 
phase-field model parameters. 

These two properties combined have made this thin- 
interface limit ideally suited to model dendritic growth in 
pure materials quantitatively at low undercooling when 
used in conjunction with efficient numerical algorithms 
[^,0. However, how to extend these results in a useful 
way to the more general case where the two phases do not 
have symmetrical properties has remained an open chal- 
lenge. Using a distinct thin-interface analysis, Almgren 
showed that the two-sided model with unequal thermal 
conductivities maps onto a modified set of sharp-interface 
equations that is plagued by finite interface thickness ef- 
fects Q. These include (i) a temperature jump across 
the interface, (ii) an interface stretching correction to the 
heat conservation condition at the interface (Stefan con- 
dition) also found previously in and (iii) a surface 
diffusion correction to the same condition. 

The same effects plague the thin-interface limit of ex- 



isting phase-field models of alloy solidification [g[|8| and 
make them inadequate to model quantitatively the low 
growth rate regime of experimental relevance. In alloys, 
(i) translates into a chemical potential jump across the 
interface associated with the well-known effect of solute 
trapping | p^ , and (ii) and (iii) modify the mass con- 
servation condition at the interface. Here, I present a 
phase-field formulation of alloy solidification that makes 
it possible to eliminate simultaneously all three effects. 
Furthermore, I demonstrate that it yields the same com- 
putational benefits as the thin-interface limit of the sym- 
metric model j^] for dendritic growth. 

For clarity of exposition, I first discuss the thin- 
interface limit of a simpler model that describes an ideal- 
ized binary alloy with parallel liquidus and solidus slopes, 
and which reduces to the standard Hele-Shaw flow prob- 
lem ||ll|] in its Laplace limit. I then consider a realistic 
dilute alloy model with unequal slopes that reduces to the 
former model in the limit where the partition coefficient 
fc ^ 1. The equations of the first model are 



Tdt(j) = W^V^<j) - [/(0) + A 5(0) u] 
dtc + V-j = 0, 



(1) 
(2) 



where c is the concentration defined as the mole fraction 
of B in a binary alloy of A and B, 

J = -AcoDq{<j))Vu - atWAcodt<j)VH\^(l)\, (3) 
and 



u = c/Acq + h{4,)/2 - (c,o + Qo)/(2Aco) 



(4) 



is a dimensionless measure of the departure of the chem- 
ical potential from its equilibrium value with u = 
in equilibrium. Cjo (c;o) are the equilibrium concen- 
trations in the solid (liquid) at a fixed temperature Tq, 
Acq = c/o — Cso, /(0) = —0^/2 -f 0'*/4 is a double well 
function with the minima </) = ±1 corresponding to the 
solid (+1) and liquid (—1), g{(t>) is an odd function of 
with 5(±1) — ±1 and vanishing first and second deriva- 
tives .g'(±l) ~ g"{±.l) = 0, and h{(j)) is an odd function 
of <j) with h{±l) = ±1 that can be chosen independently 
of g{(f>) for non- variational dynamics j^]. 
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The above model reduces to the symmetric model for 
the choice q{(j)) = 1 and at = 0. In this case, the thin- 
interface Umit of this model maps onto the Stefan prob- 
lem defined by: dtu = DV^u in both phases, the Stefan 
condition V = —D{dnu\'^ — 9„w|~), where V is the in- 
terface velocity and dnu\^ is the normal gradient of u on 
the solid (— ) and liquid side of the interface, and the 
velocity-dependent Gibbs-Thomson condition 

u = —doK — (3V, (5) 

where 

da = aiW/\ and /3 = ai [t/{W\) - a2W/D] . (6) 

The expressions for the coefficients ai and 02 are identical 
to those derived by KR and interface kinetics can be 
eliminated (/3 = 0) by choosing A = DT/{a2W^). 

Next, for the alloy case, must now be chosen to 
vary from q{—l) = 1 in the liquid to (?(+l) — DsoUd/D in 
the solid. I consider explicitly the one-sided limit where 
D solid I D 0, but the results also extend to the more 
realistic case where DsoHd/D <C 1. The essential new 
term that yields the desired thin interface limit is the 
anti-trapping mass current that corresponds to the sec- 
ond term on the right-hand-side of equation (||) , and is 
only non- vanishing in the diffuse interface region. It pro- 
duces a solute flux from the solid to the liquid along the 
direction normal to the interface that counter-balances 
the trapping current associated with the jump of chemi- 
cal potential across the interface: Au = m+ — u^, where 
correspond to the hquid and solid (— ) sides of 
the interface. Thus the anti-trapping current makes it 
possible to eliminate this jump while still leaving enough 
freedom to choose the other functions in the model to 
eliminate the corrections to the mass conservation condi- 
tion. Repeating the analyses of KR and Almgren , 
I obtain that Au vanishes if F'^ = , where 

F± - / df^[p{Mri))-p{Tl)], (7) 

and 

p(^o)= (M0o)-l + at\/2(l-0g))/g(0o). (8) 

In the above definitions, i/'oC?/) = — tanh(ry/\/2) is the 
equilibrium phase-field profile, where 77 is a coordinate 
that runs normal to the interface scaled by W , and I 
have used the identity 9,,0o = ^(1 — 0o)/V2- Next, the 
mass conservation condition has the form 

V ^ -Ddnu\+ - ciWkV - C2WDd1u (9) 

where Ci = iJ+ — , C2 = — , and 

f^°^ dTj[h{Mv))-h{Tl)], (10) 

Q±= / dr,[q{Mv))-qiTl)]- (H) 

JO 



The second and third term on the right-hand-side of 
equation (^) represent the solute redistribution due to 
stretching the interface and by diffusion along its ar- 
clength s, respectively. These two terms appear, equiva- 
lently, at two successive orders in the thin interface limit 
considered by KR, or both at second order in the distinct 
thin- interface limit considered by Almgren 

In summary, one is left with three conditions to satisfy: 
(i) F^ = F~ (chemical potential jump), (ii) iJ+ = H~ 
(stretching), and (iii) = Q~ (surface diffusion). It 
is actually possible to satisfy two of these three con- 
ditions simultaneously within the standard phase-field 
formulation without the anti-trapping current. For ex- 
ample, with at = in equation (^, the choice /i(0) — 
1 - (1 - ^)72 and g(0) = (1 - (j))/2 satisfy (i) and (iii), 
but not (ii), and it is also possible to satisfy (i) and (ii) 
but not (iii). However, test simulations show that with ei- 
ther ci or C2 non- vanishing, the morphological stability of 
the interface that sets the initial scale of the pattern be- 
comes significantly altered for computationally tractable 
choices of W ||l^ , which is also easy to check analytically 
by repeating the standard MuUins-Sekerka analysis with 
the modified boundary condition (^) . Similarly, it can be 
shown that solvability theory predictions for the dendrite 
tip become modified. Consequently, all three conditions 
must be satisfied to lift the restriction on W. 

With the anti-trapping current present, one is now free 
to make the simplest choices h{(j)) — (p and q{(j>) = (1 — 
4>)/2 that satisfy = H~ and Q'^ = Q~ , respectively. 
By choosing at = l/(2\/2), one can then reduce the func- 
tion p((/)o) to the simple formp(0o) = 0o~l that satisfies 
F^ — F^ (and a non-vanishing amount of trapping can 
also be obtained by varying a). Remarkably, with these 
choices, the thin interface limit of the one-sided model 
produces a velocity-dependent Gibbs-Thomson condition 
that is identical to the one of the symmetric model. 
Consequently the expressions for ai and 02 that deter- 
mine c?o find (3 are the same as in Ref. B: ai = I/J 
and a2 ^ {K + JF)/{2I) where / = ^^??(a^(/)o)^ 
J = .g(-l) - g{+l), F = F^ = V2\n2, and 

d7jd,,Mri)9'{Mv)) / dCMO, (12) 
-00 Jo 

where I have defined g'{(j)) = dtj,g{(t)). It also follows that 
the standard Hele-Shaw flow problem with an infinite 
viscosity contrast can be simulated by taking the 
limit 9fU — > of the present phase-field model. 

Consider now the standard one-sided dilute alloy 
model defined by the set of equations 

dtc = DV^c, (13) 
ci{l^k)V = -Ddnc\ + , (14) 
Q/c° = l-(l-fc)doK (15) 

where do = l[L\ra\{\ — k)c^] is the chemical capillary 
length, Tm is the melting temperature, L is the latent 
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heat of melting, m is the hquidus slope, k = Cg/c/ is the 
partition coefficient where q (cg) is the concentration on 
the liquid (solid) side of the interface, and solidification 
is again assumed to take place isothermally. 

To construct a thin interface limit that maps onto the 
free-boundary problem ([l^)-(|l5|), I follow the same pro- 
cedure of adding a local anti-trapping current in order to 
eliminate the jump of chemical potential together with 
the other terms. The equations are 



Tdt 



X 



(16) 



with the same continuity relation (H) as before, 



J = -Dcq{(l))Vu - atWc'^il - fc) e" V(/)/|V(/)|, (17) 



In 



2c/c° 



1 



(18) 



and where g{4>), and q{(t>) obey the same limits 

at = ±1. It is again a dimensionless measure of the 
departure of the chemical potential /i from equilibrium 
(namely here u = vo(fi — (Ie) / (RTq) where R is the rare 
gas constant and vq is the molar volume assumed to be 
constant). The logarithmic dependence of u on c is re- 
lated to the entropy of mixing in the free-energy density 
as in previous models (2[|^. The main difference here 
is again the addition of the anti-trapping mass current. 
Also, the present thin-intcrfacc analysis differs from the 
analysis of Kim et al. § that docs not consider interface 
stretching and surface diffusion, and assumes that q{(j)) 
is constant in the interface region. 

The condition for eliminating surface diffusion becomes 
now Z+ ^ Z-, where Z± = drj [M{r]) - M{±oo)] 

where M{ri) = q{4>Q{ii)) CQ{rf), and 



co(??) [l + k~{l~k)h{MT^))] /2 



(19) 



is the equilibrium concentration profile across a station- 
ary interface. Therefore choosing 
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l + k-{l- k)h{(j)) 



(20) 



satifies this condition. The choices ft.((/)) — cj) and 
at — l/(2\/2) then make the jump of u and the in- 
terface stretching term vanish up to finite corrections 
that turn out to be small in computations and will be 
discussed elsewhere. Consequently, the conditions for c 
on the two sides of the interface have the desired form, 
q/c" = 1 — (1 — k)doK — (1 — k)(3V and Cs = kci where 
the expressions for do, (3, oi, and 02 are again identi- 
cal to those for the symmetric model quoted earlier here. 
Therefore, p can again be made to vanish. The dilute al- 
loy model is easily shown to reduce to the parallel slope 
model in the limit fc — > 1 by making the change of variable 
U = (e" — 1)/(1 — fc). A small solid diffusivity can also be 



modeled by adding (1 + cj)) D soUd / i'^D) to the right-hand- 
side of equation (po|), which has a negligible effect on the 
thin-interface limit in typical alloys where solid diffusion 
is substitutional (D^oiid/D^ IQ-^) 

Finally, the model is straightforward to extend to di- 
rectional solidification with the standard frozen temper- 
ature approximation T{z) = Tq + G{z — Vpt) that yields 
(for P — 0) the interface condition 

c//c? = 1 - (1 - k)doK - (1 - fc)(z - Vpt)/lT, (21) 

where Vp is the pulling velocity of the sample along the 
z-axis. It = jn^Kl — k)c^ /G is the thermal length , G is 
the temperature gradient, and c" — Coo/k where Coo = 
c{z = -l-oo). It simply suffices to substitute e" by -I- 
(1 — k){z — Vpt)/lT in equation (p^). 

The convergence of the model was examined by carry- 
ing out two-dimensional simulations of isothermal den- 
dritic growth. These simulations are directly analogous 
to the ones carried out previously to test the thin in- 
terface limit of the symmetric model Crystalline 
anisotropy was included by generalizing equation ( p^ ) to 
a standard anisotropic form 



-(9) dt(t> = 



f{<t>) 



X 



■5(0)(e--l) 



(22) 



+ V • {W{ef^(j)) - (W{Q)W'(Q)dy(l)) 
+ dy{W{Q)W'{Q)d^^), 

where <d = a,icta,n{dy(j)/dx4>) is the angle between the 
direction normal to the phase-field interface and the 
X-axis. As a result, the anisotropic form of do and 
(3 become do{e) = ai[Wie) + W"{e)]/X and (3{e) = 
ai [T{e)/{W{9)X) - a2W{9)/D] in the interface condition 
where ai = 0.8839, and 02 = 0.6267 for the common 
choice d^g{<f>) = {l — cjP)'^. Furthermore, I chose the stan- 
dard form of four-fold anisotropy W{9) = Was{0) with 
as{0) ~ I + £4 cos 46* and made (3{6) vanish by letting 
T{e) = rasiB)^ and A = DT/{a2W'^). 

I simulated equations (||) and ( p2| ) with j and u de- 
fined by equations (|l^) and (|l^), respectively. I com- 
pare the results of the present model with at = l/(2\/2) 
and q{(j)) defined by equation (pO|), to the more stan- 
dard choice (i.e. similar to previous models) that has 
no anti-trapping current (at = 0) and uses the simplest 
scaled diffusivity function (/((/)) = (1 — (/))/2; h{(j)) = (j) and 
d^g{(t)) = (1 — (/>^)^ in both models. The former model has 
the desired thin-interface limit with local equilibrium at 
the interface, whereas the latter has both a chemical po- 
tential jump and surface diffusion. I used a simple finite- 
difference Euler method with Ax — 0.4 and Ai = 0.008, 
W = T = 1, £4 = 0.02, k = 0.15, and the scaled supersat- 
uration ft — (c^ ~ Coo)/[c°(l — k)] = 0.55 where Cqo is the 
initial alloy concentration. In all simulations, the initial 
condition consisted of a circular seed of radius r — 22 da, 
u — ln(l — (1 — k)Vl), and c defined by equation ( p^ that 
varies from Coo in liquid to fccoo in solid. 
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FIG. 1. Plots of scaled dendrite tip velocity Vdo/D vs 
scaled time tD/dg for at = l/{2\/2) and q{(j)) given by Eq. po| 
(present), and at — and q{(l}) = (1 — 0)/2 (standard). 
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FIG. 2. Plots of solute profiles in the solid along the cen- 
tral dendrite axis, and comparison with the Gibbs-Thomson 
relation for the present model (short-dash line). 



The dimensionless dendrite tip velocity Vdo/D is plot- 
ted vs the dimensionless time tD/d^ for the two models 
and two different ratios of do/W in Fig. |l|. Since the 
simulation time scales ~ (do/W)^, the runs with da/W 
twice smaller are ~ 32 times shorter. Furthermore, I 
have plotted in Fig. |^ the scaled concentration Cs{x)/c'i 
in the solid vs the scaled position x/do along the cen- 
tral dendrite axis for the two different models. Plots in 
Figs. |l|-^ must superimpose when results are converged. 
These plots show, as expected, that the present model is 
well converged in this range of do/W that is comparable 
to the one studied in the symmetric model whereas 
the standard model is not. This is especially true for 
the microsegregation profile that is still far from being 
converged in the latter model, even for the largest ratio 
do/W = 0.544. In contrast, this profile is already well 
converged for a twice smaller ratio in the present model; 
it agrees, self-consistently, within a few percent with the 
Gibbs-Thomson relation Cs{x)/cf = k[l — (1 — k)do/p] 
where p is the dendrite tip radius in the simulation. It 
will be shown elsewhere that this dramatic difference of 
convergence for microsegregation is due to the fact that 
the amount of solute trapped ^ plnp for small velocity 
[p = WV/D <C 1) in the standard model. 

The present results demonstrate that the phase-field 
method can be successfully extended to model quantita- 
tively microstructural pattern formation in alloys with a 
realistic solid diffusivity. For this important application, 
it is potentially more advantageous than the level set 
method p3| since it does not require the explicit compu- 
tation of the interface velocity. These results also revive 
the hope to extend the phase-field method to model ac- 
curately a wide range of other interfacial patterns with a 
strong asymmetry between phases. 
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